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An approach is proposed to improve the efficiency of fourth-order algorithms for numerical in- 
tegration of the equations of motion in molecular dynamics simulations. The approach is based 
on an extension of the decomposition scheme by introducing extra evolution subpropagators. The 
extended set of parameters of the integration is then determined by reducing the norm of truncation 
terms to a minimum. In such a way, we derive new explicit symplectic Forest-Ruth- and Suzuki-like 
integrators and present them in time-reversible velocity and position forms. It is proven that these 
optimized integrators lead to the best accuracy in the calculations at the same computational cost 
among all possible algorithms of the fourth order from a given decomposition class. It is shown 
also that the Forest-Ruth-like algorithms, which are based on direct decomposition of exponential 
propagators, provide better optimization than their Suzuki-like counterparts which represent com- 
positions of second-order schemes. In particular, using our optimized Forest-Ruth-like algorithms 
allows us to increase the efficiency of the computations more than in ten times with respect to 
that of the original integrator by Forest and Ruth, and approximately in five times with respect to 
Suzuki's approach. The theoretical predictions are confirmed in molecular dynamics simulations of a 
Lennard- Jones fluid. A special case of the optimization of the proposed Forest-Ruth-like algorithms 
to celestial mechanics simulations is considered as well. 
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I. INTRODUCTION 



Modelling various physical and chemical processes in 
molecular dynamics (MD) simulations wc come to the ne- 
cessity to integrate the equations of motion for a many- 
body system of interacting particles. A lot of numeri- 
cal algorithms have been devised and implemented over 
the years to perform such an integration. The tradi- 
tional high-order explicit Runge-Kutta (RK) and implicit 
predictor-corrector (PC) schemes were applied in 
early investigations. Is was soon realized that the ex- 
tra orders obtained in these schemes are not relevant 
since the truncation errors accumulate drastically on MD 
scales of time |^ . This high instability restricts the appli- 
cation of RK and PC integrators in long-term MD simu- 
lations to very small time steps only, and, thus, reduces 
significantly the efficiency of the computations. In addi- 
tion, the RK and PC algorithms produce solutions which, 
unlike exact phase trajectories, are neither symplectic nor 
time reversible. 

In 1990, a new approach to the integration of motion 
in many-body systems has been proposed Within 
this approach, the time propagation is carried out on the 
basis of exponential decompositions of evolution propaga- 
tors. The main advantage of the decomposition method 
is that for an arbitrary order in the time step it allows to 
construct algorithms which are exactly symplectic and 



time-reversible. The preservation of symplecticity and 
reversibility appears to be very important because, as is 
now well established, this closely relates to the stability 
of an algorithm ||^ . Another nice property of the decom- 
position integration is its explicitness and simplicity in 
implementation. This is in a sharp contrast to implicit 
time-reversible symplectic algorithms obtained recently 
pO| within the RK approach, where cumbersome sys- 
tems of coupled nonlinear equations must be solved by 
iteration at each step of the integration process. 

Nowadays, the decomposition method should be con- 
sidered as the main tool for construction of efficient inte- 
grators of motion in classical as well as quantum systems 
Q. Modified versions of this method have also been in- 
troduced. In particular, it was shown that for atomic 
systems with long-range interactions the efficiency of the 
integration can be improved by additionally splitting the 
Liouville operator into slow and fast components |ll^,f2| . 
In such a multiple scale propagation, the slow subdynam- 
ics is treated in a specific way using larger step sizes in 
view of the weakness of long-range forces. The fasted 
motion, caused by the interactions at short interparti- 
cle distances, remains to be integrated with the help of 
original decomposition algorithms. The question of how 
to derive higher-order integrators by composing lower- 
order decomposition schemes has been considered as well 
13|Jl4|. Moreover, it has been shown recently how to 
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adapt the decomposition approach to integrate not only 
translational motion in atomic Uquids, but also simulate 
more complicated molecular and spin systems at the pres- 
ence of orientational degrees of freedom |]l5|-|l7| . 

The main attention in previous studies has been di- 
rected to obtain algorithms which require a minimal 
number of force evaluations per time step. For in- 
stance, the well-known velocity- and position- Verlet in- 
tegrators |ll|,|l^ correspond to a second-order decompo- 
sition scheme with one force evaluation per step. The 
fourth-order algorithm by Forest and Ruth Q presents a 
scheme with three such force recalculations. Sixth-order 
schemes can be reproduced beginning from seven 
evaluations of force for each particle during each time 
step. It is worth mentioning that sixth- and higher-order 
schemes involve too large number of force recalculations 
and generally are not recommended to be used in MD 
simulations. The reasons are that for a system with a 
great number of particles, the force evaluations consti- 
tute the most time-consuming part of the computations 
and that the truncations errors decrease with increas- 
ing the order of the scheme much slower than the total 
number of these evaluations. Such high-order algorithms 
may be efficient only for systems composed of a few bod- 
ies when a very high precision in determination of phase 
trajectories is desirable, for instance, in astronomical ap- 
plications. In the case of MD simulations, the preference 
should be given to more simple second- or fourth-order 
schemes. For many MD applications to achieve a required 
level of accuracy, the fourth-order integrators appear to 
be more efficient than second-order schemes. 

The minimal numbers of force evaluations per time 
step do not guarantee, however, the optimization with 
respect to the overall number of such evaluations which 
are necessary to perform during a fixed observation time 
interval. This is so because schemes, in which these num- 
bers exceed the minimum, can be used with larger sizes 
of the time step in order to obtain the same accuracy in 
solutions. Only very few papers were devoted 

to derivation of such extended schemes. In particular, 
Suzuki et al. have pointed out that the fourth-order al- 
gorithm by Forest and Ruth may not lead to optimal 
performance since during the evaluation it involves time 
coefficients which are larger in magnitude than the size 
of the initial time step Q. As a result, the original algo- 
rithm has been replaced by a new fourth-order integra- 
tor composing five, instead of three, second-order Verlet 
schemes. However, the question of how this replacement 
influences on the efficiency of the computation has not 
been studied. Extended Suzuki-like schemes were also 
the subject of investigations by Kahan and Li |l9j. Com- 
posing second-order algorithms, they introduced higher- 
order integrators with time coefficients chosen to provide 
a minimum for the maximal one among them in magni- 
tude or for the sum of absolute values of these coefhcients. 
Again, it has been unknown whether the increased num- 
ber of force evaluations in their integrators is compen- 
sated by the possibility of using larger step sizes or not. 



In Ref. [gOj] new explicit velocity- and position- Verlet-like 
algorithms of the second order were proposed to integrate 
the equations of motion in many-body systems. These 
algorithms were derived on the basis of an extended de- 
composition scheme at the presence of a free parameter. 
The nonzero value for this parameter was obtained by 
reducing the influence of truncated terms to a minimum. 
As a result, the new optimized algorithms appear to be 
more efficient than the original Verlet versions which cor- 
respond to a particular case when the introduced param- 
eter is equal to zero. 

In the present paper we propose a consequent approach 
for improving the efficiency of fourth-order decomposi- 
tion integration. As a result, we derive new Forest-Ruth- 
and Suzuki-like integrators by explicitly reducing the in- 
fluence of truncation terms to a minimum. Although this 
requires (as in any extended scheme) extra force evalua- 
tions per time step, the overall number of force recalcula- 
tions keeps unchanged due to the use of larger step sizes. 
At the same time, the resulting precision increases sig- 
nificantly with respect to that of the original algorithms 
by Forest and Ruth, and Suzuki. It is demonstrated also 
that previous criteria used for optimization of decompo- 
sition algorithms may have nothing to do with providing 
the best performance of the computations. 



II. DECOMPOSITION INTEGRATION 

We will deal with a classical iV-body system described 
by the Hamiltonian 

1=1 i^j 

where and Vj denote the position and velocity, respec- 
tively, of particle i with mass m, and firij) is the inter- 
particle potential of interaction with = jr^ — Vj\. The 
equations of motion for such a system can be cast in the 
following compact form 

^ = [poH]^Lpit). (2) 

Here p = {r^, v^} is the full set (i = 1, 2, . . . , iV) of phase 
variables, [o] represents the Poisson bracket and 

2 = 1 

is the Liouville operator with fi = — X^^^yj) (p'{rij)rij /rij 
being the force acting on particles due to the interaction. 

If an initial configuration p(0) is provided, the unique 
solution to Eq. (2) can be presented as 

p{h) = e^'^piO) = e(^+^)'V(0) , (4) 
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where h denotes the time step, and the Liouville op- 
erator L = A + B has been spUt into the free-motion 
A — v-d/dr and potential B — i/m-d/dv parts with 
V = {vi}, r = {r^}, and f = {fi}. The significance 
of such a spUtting will be understood below. Then the 
evolution of the system can be investigated during arbi- 
trary times t by repeating the single-step propagation, 

p{t) = (e^'') p(0) = (e(^+-»"')'p(0), where I = t/h is 
the total number of steps. 

Of course, solution (4) is quite formal because the ex- 
ponential propagator e^'' does not allow to be evaluated 
exactly at any h (solutions in quadratures are possible 
only for N = 2 that is not relevant for our case of many- 
body systems when A'' ^ 1). However, at small enough 
values of h, the total propagator can be decomposed 
using the formula 



,(A+B)/i+0(/i^+i) 



n 



ph^Bbph 



(5) 



The decomposition method is quite general to build 
numerical integrators of arbitrary orders. In particular, 
the second-order {K — 2) velocity- Verlet (VV) algorithm 
[|ll|,|8[ is obtained from Eq. (5) at P = 2 and oi = 0, &i = 
62 = 1/2, aa = 1, i.e., e^A+B)h+o(h^) ^ ^Bh/2^Ah^Bh/2^ 
The case when the operators A and B are replaced by 
each other {A ^ B) is also possible, and we come 11| to 
the position- Verlet (PV) integrator 



JA+B)h+0(h^) 



(7) 



corresponding to the choice ai = a2 — 1/2, bi — I, and 
&2 = 0. The fourth-order {K — 4) algorithm by For- 
est and Ruth (FR) is immediately reproduced from 
Eq. (4) at P = 4 



g(A-hB)/i+0(h=) ^ ^A0^^Beh^A{l-e)^^B{l-2e)h ^ 

^A{i-e)^^B9h^Ae^ 



(8) 



The coefficients ap and bp in this formula should be cho- 
sen in such a way to provide the highest possible value 
for > 1 at a given integer number P > 1. Then in- 
tegration (4) can be performed approximately with the 
help of Eq. (5) by neglecting truncation terms 0{h^^^). 

The main advantage of the above decomposition is that 
the exponential subpropagators e^"^ and c^'^, appearing 
in the right-hand-side of Eq. (5), are analytically inte- 
grable. Indeed, as can be shown readily, 



e-^V EE e^'^/*^{r, v} = {r + rv, v}, 
e^^p = ef/"-^/^^^{r, v} = {r, v + rf /m} 



(6) 



with T being equal to aph or bph. That is why we ex- 
plicitly presented L as the sum of free-motion A and 
potential B terms. Another important feature of the 
decomposition integration is that it leads to symplec- 
tic trajectories. This is so since Eq. (6) represents, in 
fact, simple shifts of position and velocity, and these 
shifts do not change the volume in phase space. The 
time reversibility S{—t)p{t) — p{0) of solutions (follow- 
ing from the property S~^{t) = S{—t) of evolution op- 
erator S{t) = e^*) can also be reproduced by impos- 
ing additional constraints on the coefficients ap and bp, 
namely, oi — 0, Cp+i = ap_p+i, and bp — 6p_p+i, or 
Up = ap-p+i, and bp — ap^p a.t bp — 0. Then the sub- 
propagators e"^"^ and e^"^ will enter symmetrically in the 
decompositions and, thus, provide automatically the re- 
quired reversibility. Note also that such constraints lead 
to automatic disappearing even-order terms oc h'^^ in the 
function 0(/i^+^) for any fc > 0. For this reason, the 
order K of time-reversible algorithms may accept only 
even numbers, (K = 2,4,6,...). The cancellation of 
odd-order terms oc h'^'^^^ in 0{h^'^^) up to a required 
finite number k will be provided by fulfilling a set of ba- 
sic conditions for Uk and bk- For example, the condition 
Sr=i '^p ~ S^=i = 1 is necessary to cancel the first- 
order truncation uncertainties. 



with ai = 04 = e/2, a2 = 03 = (1 - e)/2, 61 = 63 = 0, 
62 = (1 - 20), 64 = 0, and 6 = 1/(2 - ^/2) « 1.3512. 
Propagation (8) can be related to the position version of 
the FR integrator (PER) , because putting formally — 
it transforms to the second-order PV algorithm (7). For 
A ^ B, Eq. (8) will represent the velocity FR counter- 
part (VFR), which can be derived directly from Eq. (5) 
at ai = 0, 02 = 04 = 0, aa = (1 — 20), 61 = 64 = 6/2 and 
&2 - 63 = (1 - 0)/2. 



III. OPTIMIZATION OF FOURTH-ORDER 
ALGORITHMS WITHIN DIRECT 
DECOMPOSITION 



Let us consider now an extended decomposition scheme 
of the fourth order by allowing to accept a value for P 
which exceeds the necessary minimum (P = 4) on unity, 
i.e., letting P = 5. Remember that we cannot choose 
the number P to be too big, because this results in too 
larger number, namely P — 1, of expensive force recalcu- 
lations. Chossing P = 5 we hope simply to reduce the 
truncation errors 0{h^) significantly in a little additional 
computation cost, rather than to increase the order of the 
decomposition scheme (note that sixth-order integrators 
are derivable [Q beginning up from P = 8). 

For P = 5, the extended decomposition can be pre- 
sented in the form 



g(A+i3)/i+C3/rVC5/i"+0(/i^) ^ gi3ChgA(l-2A)f gSx'ix 
AXh^B(l-2{x+0)hQA\h^Bxh^A(l-2X)^^B^h 



(9) 



e e 



'2e" 



following from Eq. (5) at ai = 0, 61 = 65 = ^, 02 — 
as = (1 - 2A)/2, 62 = 64 = X: 03 = 04 = A, and 
63 = 1 — 2(x + 0- Here the symmetry of time coeffi- 
cients and the condition J2p=i = J2p=i bp — I have 
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already been taken into account. Again, the propagation 
with A <-> i? is also acceptable, and then ai = = ^, 
61 = 65 = (1 - 2A)/2, 02 = 04 = X, h = bi = A, 
03 = 1 — 2{x + £,), and 65 = 0. The operator C3, appear- 
ing in the left-hand-side of Eq. (9), is responsible for the 
cancellation of third-order truncation uncertainties. The 
explicit expression for it is 

- A, x) [A, [A, B]] + X,x)[B, [A, B]] , (10) 

where [ , ] denotes the commutator of two operators, and 



A, x) = + Ax(l - X - 2^ + I - ^ 



(11) 



So that at C3 = 0, i.e. when q;(^, A, x) = and 
/3(^, A,x) = 0, formula (9) represents a whole family of 
symplectic time-reversible integrators of the fourth-order. 

A particular member of the above family can be ob- 
tained by choosing corresponding values for ^, A, and x- 
As far as there are three parameters and only two con- 
straints, a = and /3 = 0, one from these parameters, ^ 
say, can be treated to be free. Then, for example, putting 
5 = 0, Eq. (9) reduces to the original FR algorithm (8) in 
position or (when A ^ B) velocity forms. The extended 
(when 5 7^ 0) propagation will require already four, in- 
stead of three, force recalculation per time step. How- 
ever, having a room in varying ^, we can overcompensate 
the increased computational efforts by minimizing the 
fifth-order truncation uncertainties C^h^. 

In order to show this, let us analyze in detail the in- 
fluence of these uncertainties on the result. Expanding 
both the sides of Eq. (9) into Taylor's series with respect 
to h, one finds 

C, = 71 [A, [A, [A, [A, B]]]] +72 [A, [A, [B, [A, S]]]] + 

73 [i?, [A, [A, [A, B]]]]+^4B, [B, [B, [A, B]]]]+ (12) 
75 [B, [A, [A,B]]]]+^e[A, [B, [B, [A,B]]]] , 

where explicit expressions for 7-multipliers are: 



71 



72 



73 



74 



75 
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Assuming that all the fifth-order commutators in ( [T^ ) 
are nonzero valued, the norm of C5 with respect to fifth- 
order commutators arising in Eq. (12) can be written as 
follows 



7(C, A, x) - ^ 7? + 7| + 7| + 7| + 7| + 76 • (13) 

Then the norm of local uncertainties C^ph^ appearing in 
phase trajectory p during a single-step propagation given 
by Eqs. (4) and (9) can be expressed in terms of 7 and 
h as g — jh^. During a whole integration over a fixed 
time interval t, the total number I of such single steps is 
proportional to h~^. As a result, the local fifth-order un- 
certainties will accumulate step by step leading a.t h 
to the fourth-order global errors F = gh~^, i.e.. 



F(/i,e,A,x)=7(e,A,x)/i' 



(14) 



Extended propagation (9) can now be optimized with 
respect to time coefficients A, and x by finding the 
global minimum for the function 7(^1 A, x), provided 
a — and /3 = 0, i.e., solving the system of equations 



«(C,A,x) =0, 
/3(5,A,x) =0, 
7(^7 A, x) = min (global) . 



(15) 



A way to simplify the problem is to solve analytically the 
first two equations of (15) with respect to x and 5, 



4A - 8A2 ± J2X(-1 + 8A - 24A2 + 24A3) 
X(A) = 



12A - 96A3 + 96A4 



(16) 



e(A) = --4A2x, 
6 



then substitute expressions (16) into the third equation, 
and find numerically the minimum considering already 
7(5(A), A, x(A)) as a function of only one variable A. The 
result within sixteenth significant digits is 



C = +0.1720865590295143E-h00 
A = -0.9156203075515678E-01 
X = -0.1616217622107222E-h00. 



(17) 



The global minimum of 7(5, A, x), corresponding to 
solutions (17) found (they satisfy Eq. (16) at sign 
in the right hand side of expression for x(A)), consists 
7mhf'" ~ 0.00092, where the superscript EFRL refers to 
the extended FR-like integration (9). On the other hand, 
the value of 7 corresponding to usual FR scheme (8), i.e. 
when 5 = 0, A = (1 - 6')/2 and X = 6* (then Eq. (9) 
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reduces to Eq. (8)), is equal to 7fr ~ 0.039. We see, 
therefore, that applying the extended integration allows 
one to decrease the truncation errors approximately in 
iFR/Trnfif"^ ~ 42 times. Taking into account that such 
a decreasing has been achieved increasing the number 
of force evaluations per step from three to four, the ex- 
tended propagation must be performed with step sizes 
which are in factor 4/3 higher than those of the FR algo- 
rithm, in order to provide the same number of total force 
recalculations during the fixed overall interval of integra- 
tion. Thus, we will come to more efficient calculations 
if the following inequality rJ^J,fL(4/),/3) < Tpnih) takes 
place. In view of Eq. (14), such an inequality can be 
rewritten as 7FR/7m™^ > (4/3)'' « 3.16, so that it is 
fulfilled completely in the optimization regime. In par- 
ticular, 



TFR{h) 



0.075 



(18) 



indicating that the global errors can be reduced more 
than in 10 times with respect to the FR integration with- 
out spending any additional overall computational costs. 

The proposed procedure can be used, in principle, for 
the decompositions of arbitrary two noncommutativc op- 
erators A and B to achieve the best performance. Note 
that the necessity in these decompositions may arise not 
only when considering the integration of motion in classi- 
cal systems, but also in mathematical and quantum me- 
chanical calculations. In some cases, further improve- 
ment of the efficiency of the computations may also be 
possible. For instance, an extra optimization of the de- 
composition scheme can be achieved in celestial mechan- 
ics applications due to a specific character of motion in 
the solar system (see Appendix). 

In the case of MD simulations with vclocdty-indepcn- 
dent forces, an additional optimization can be carried 
out as well using specific properties of operators A and 
B. Taking into account explicit expressions for these 
operators, it can be verified readily that two of the 
sixth five-order commutators vanish in Eq. (12), namely, 
[B, [B, [B, [A, B]]]] = and [A, [B, [B, [A, B]]]] = 0, when 
time propagation is performed with the help of extended 
scheme (9). This scheme, by analogy to Verlet integra- 
tors, we will refer to the velocity version of the EFRL 
algorithm and abbreviate in our notations as VEFRL. 
Then, letting formally 74 = and 76 = to cxchide 
the above zeroth commutators in Eq. (13), and resolving 
problem (15) on condition minimum yields 



^ = -F0.1644986515575760E-K00 
A = -0.2094333910398989E-01 
X = -F0.1235692651138917E-I-01 . 



(19) 



This leads to the global minimum 7^?f ~ 0.00065 

(which is achieved at sign "— " in Eq. (16)), while the 
value of 7 corresponding to usual FR integration (8) 
reduces (at 74 = and 76 = 0) to 7fr ~ 0.028. 



Again, the truncation uncertainties decrease approxi- 
mately in the same factor, iFR/Tmkf^'" ~ 43, and 
rlfJ^^{4h/3)/TMh) « 0.073. 

The position EFRL integration (PEFRL) is obtained 
from Eq. (9) by replacing A ^ B. In such a case, the 
zeroth five-order commutators are [A, [A, [A, [A, B]]]] = 
and [B, [A, [A, [A, B]]]]. They can be excluded in Eq. (13) 
by putting 71 = and 73 0. Then solutions to system 
(15) transform into 



^ = -h0.1786178958448091E-|-00 
A = -0.2123418310626054E-I-00 
X = -0.6626458266981849E-01 



(20) 



(with sign "+" in Eq. (16)), and the minimum is 
Imirf^^ ~ 0.00061. The value of 7 corresponding to 
original scheme (8) is now (when 71 = and 73 = 0) 
equal to 7fr ~ 0.038. so that jFRh^fJ^^ « 62 and 
r^™(4V3)/rFR(/i) « 0.051. 

In view of Eqs. (4), (6), and (9), more explicit ex- 
pressions for the single-step propagation of position and 
velocity from time t to t+h within the optimized VEFRL 
algorithm are: 

VI = v(f) + ^{[rmh 
n = r{t)+vi{l - 2X)h/2 
V2 = vi ^f[ri\xh 
r2 = ri -I- V2A/1 

V3=V2 + if[r2](l-2(x + 0)/i (21) 

r3 = i"2 -I- V3A/1 

V4 = V3 ^t[r3]xh 
r{t + h) = r3+ V4(l - 2X)h/2 
v(f + /7,)=V4 + ^f[r(t + /i)]^ft 

where the values for ^, A, and x should be taken from 
Eq. (19). The optimized PEFRL algorithm (when A ^ 
B in Eq. (9)) reads: 



ri 


= r(i) + v{t)^h 




Vl 


= v(i) + if[ri](l- 


- 2X)h/2 


r2 


= ri +vixh 






= Vl + if[r2]Xh 






= r2 + V2(l - 2(x- 


+ 0)h 


V3 


= V2 + ^i[r3]Xh 




r4 


= rs + V3X/1 




v(t + h) 


= V3 + ^f[r4](l- 


2X)h/2 


r{t + h) 


= r4 -h v(t -1- h)]^h 





(22) 



and the parameters ^, A, and x should accept their val- 
ues from Eq. (20). The algorithms arc explicit, simple 
in implementation, and require only slight modification 
with respect to the original FR scheme. 
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IV. OPTIMIZATION BY COMPOSING 
SECOND-ORDER SCHEMES 

Another way to construct high-order algorithms con- 
sists in composing lower-order schemes. In particular, 
employing second-order Verlet integrator (7) , the compo- 
sition can be performed 0| for an arbitrary higher order 
K >2 as 



Q 

n 

9=1 



9=1 



(23) 



Here the coefficients dg are chosen by providing a max- 
imum for K at a given number Q > 1. This leads 
to the necessity of fulfilling a set of order constraints, 
and, for instance, the condition dq — 1 must 

be satisfied to avoid the first-order truncation terms in 
0{h^^^). For time-reversible compositions, the coeffi- 
cients dq should appear symmetrically, i.e., satisfy the 
property dq = dg-q+i. Then, as in the case of direct de- 
composition (5), even-order terms will be absent in the 
function 0{h^'^^), and the order of reversible composi- 
tion schemes will also accept only even numbers. 

Fourth-order {K = 4) composition integrators are 
derivable from Eq. (23) at Q > 3. For example, putting 
Q = 3 as well as di = d-s = and d2 — 1 — 29 one comes 
to the scheme 

eiA+B}h+o(h'') ^ S2{0h)S2{{l - 2e)h)S2{0h) (24) 

which coincides entirely with the FR algorithm (as can 
be seen by comparing with Eq. (8)). Thus, in this par- 
ticular situation, the direct decomposition and second- 
order-based composition approaches exhibit to be iden- 
tical (but with increasing K oi Q both the approaches 
will lead to different results). We will now consider the 
construction of extended (when Q > 3) composition al- 
gorithms of the fourth order. 

It can be shown that no real solutions exist for dq at 
(5 = 4 and K = A. So putting Q = 5, one obtains the 
simplest extended fourth-order composition formula 



^{A+B)h+C3h^+Ci,h-' + 0{h'') ^ 

S2{Ch)S2iXh)S2iil - 2(e + X)h)S2iXh)S2m 



(25) 



following from Eq. (23) at c?i = ds = S,, d2 = d^ — A, and 
d3 = l-2{^ + X), where 

C3 = -a(e, A) [A, B]] + ^[B, [A, B]]) (26) 



with 



,{^,X)=2e + 2X' + {l-2{^ + X)y 



Formula (25) constitutes a family of composition time- 
reversible integrators of the fourth-order, provided C3 = 



0, i.e. a(^. A) = 0. Any one from the two param- 
eters ^ and A can be chosen, in principle, arbitrarily 
since we have only one constraint. For ^ = and 
X ~ 0, we come to usual FR integration (8). Suzuki 
[0, considering an extended scheme like (25), has im- 
posed the additional constraint ^ = A = "iJ, and obtained 
^ = 1/(4 - v/4) w 0.41449. We will show below that al- 
though Suzuki's approach leads to an increased efficiency 
with respect to the original FR scheme, it is not the best 
choice for fourth-order integration. 

The operator C5, which forms the fifth-order term of 
truncation uncertainties in composition (25), can again 
be presented in form of Eq. (12), where now 
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(28) 



with 

s(C, A) = 1 - 10A(1 - 8^ -f 24^^ - 32^3 ^g^4) _^ 
40X^(1 12^2 _ 8^3) _ 80X^(1 - 4^ + 4^^) + 

80A'*(1 - 20 - 30A^ - 10^1 - 4^ + 8^^ - 8^^ + 3^^) 
and 

wiC, A) = (-A(l - 14^ + 54^2 _ 80^3 ^ 40^4^ _^ 

A2(7 - 54^ + 120^2 _ 80^3) _ ^3(^^j _ 74^ ^ 74^2) ^ 
A*(17 - 34^) - 6A5 - e + 7C^ - 17^^ + 17^* - 6^^)76 . 

It is worth remarking that the operators C3 and C5 can- 
not be reduced to zero simultaneously at Q = 5, since we 
cannot satisfy all the three equations a = 0, s = 0, and 
w = having only two parameters ^ and A (such a re- 
duction, resulting in sixth-order composition integrators, 
is possible within real coefficients dq beginning up from 
Q = 7§). ^ 

Our task is to minimize the norm 



7(e,A) = ^ 7? + 7| + 7| + 7| + ll + ll 



(29) 



of C5 in space of fifth-order commutators, provided C3 
0, i.e., to solve the problem 



a(OA) =0, 

7(0 A) = min (global) , 



(30) 



The simplest way to do this is to transform the first of 
two equalities of system (30) from the cubic (see Eq. (27)) 
to square equation replacing the variable ^ by the new 
independent quantity x — 1~2(^-|-A), and find solu- 
tions for A considering x as a parameter. As a result, we 
obtain 



(27) A(x) = 



3(1 - x)2± V3(-l + 4x - 6x2 - 12x3 + 15^4) 



12(1 -X) 



^(x) = (l-2A(x)-x)/2. 



(31) 
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Then, substituting these expressions into the second 
equation, we come to the function 7(C(x)j -^(x)) = 7(x) 
which depends aheady on only one variable x- The 
global minimum of this function is achieved at x = 
-0.7269082885036828E + 00 (with sign "+" for A(x) in 
Eq. (31)) and consists 7^^?^ « 0.0011, where the super- 
script ESL refers to the extended Suzuki-like integration 
(25). So that according to Eq. (31), the corresponding 
solutions are 



^ = 0.3221375960817984E+00 
A = 0.5413165481700430E+00 . 



(32) 



As was mentioned above, Suzuki has used the param- 
eters ^ = A = 1? « 0.41449 in his integration. This 
corresponds to the value 73 « 0.0015 of function (29) 
that is approximately in factor 1.5 higher than its opti- 
mized counterpart 7^fn ~ 0.0011. Therefore, his choice 
was not optimal. But the main conclusion we may make 
in this context is that the original Suzuki approach as 
well its optimized ESL version are less efficient than the 
EFRL integration described in the preceding section. In- 
deed, taking into account the value 7i^i™'^ ~ 0.00092 and 
the fact that the Suzuki approach requires up five (see 
Eq. (25)), instead of four, force evaluations per time step, 
we obtain in view of Eq. (14) that efficiency of the EFRL 
integration is in factor 



mm \^ / I I] 



ESL 



^ESL 

min 

T.EFRL/L^ ~ EFRL 

mill V / / min 



(33) 



better with respect to the ESL scheme, and, thus, ap- 
proximately in 5 times higher with respect to the original 
approach by Suzuki. 

Note that results (32) and (33) correspond to a gen- 
eral case of integration of motion when forces may ex- 
plicitly depend, rigorously speaking, on velocities. For 
velocity-independent accelerations, we can exploit addi- 
tional properties of operators A and B, to improve the 
efficiency of the compositions. Mention that such proper- 
ties are [B, [B, [B, [A, B]]]] = and [A, [B, [B, [A, B]]]], or 
[A, [A, [A, [A, B]]]] - and [B, [A, [A, [A, B]]]] atA^B, 
and they can be taken into account by putting formally 
74 = and 76 = 0, or 71 = and 73 = in Eq. (29). 
Then within our basic definitions for A and B, Eq. (25) 
will lead to the position version of the ESL integration 
(PESL). For this version the optimal solutions we found 
are 



C = 0.3162227486360109E+00 
A = 0.5521563637246984E-h00 



(34) 



with the minimum value 7mkf'" ~ 0.00109 that corre- 
sponds to X = -0.7367582247214187E+00 (and sign "-h" 
in Eq. (31)). Remembering that for the position version 
of the EFRL integration was 7mf/^^ w 0.00061, one ob- 
tains rPfSL(5/j/4^/rPEFRL(^^ _ 4^ whereas with respect 

to the position version of original Suzuki's approach (PS) 



for which 7ps ~ 0.0014, the efficiency increases more 
than in 5 times. It is necessary to point out that so- 
lutions (34) correspond to a local minimum of 7 (when 
74 = and 76 — 0). The global minimum is achieved 
at C = 1.453335388755048, A = -2.154220603244978, 
with X = 2.401770428979862 and consists about 0.00096 
that is only slightly smaller than the presented above 
value 0.00109. But now all the time coefficients appear 
to be significantly larger than their counterparts (34) . In 
such a situation, the preference should be given to the lo- 
cal minimum, because it leads to a significantly smaller 
value for the norm of higher-order truncation uncertain- 
ties 0{h7), and thus the precision of composition (25) 
will be much less sensitive to increasing the size of the 
time step h (see the next section). For A^ B, Eq. (25) 
will correspond to the velocity version of the ESL inte- 
gration (VESL) and its solutions (at 71 = and 73 = 0) 
to Eq. (30) take the values 



i = 0.3226106225667342E+00 
A = 0.5404642725582767E-H00 



(35) 



with X = -0.7261497902500218E-h00 (sign "-f-" in 
Eq. (31)) and the global minimum 7^?^^^ w 0.00107. 
At the same time, 7^^^^'^ ~ 0.00065 and for the veloc- 
ity version of Suzuki's approach (VS) we obtain 7vs ~ 
0.0014. So that T^f^^{5h/4:)/r^J^^{h) « 4, and again 
the ratio of efficiencies of the VEFRL and original VS 
scheme is equal approximately to 5. 



V. APPLICATION TO MOLECULAR DYNAMICS 
SIMULATIONS 

In order to verify theoretical predictions presented in 

sections III and IV, wc have applied the velocity and po- 
sition versions of our EFRL and ESL algorithms to MD 
simulations of a Lennard- Jones fluid (LJ), and compared 
the results with those of the original FR and Suzuki's 
approaches. The system under consideration was a col- 
lection of N = 256 particles interacting through a shifted 
LJ-likc potential, (f{r) = $(r) — <I>(rc) at r < Tc and 
ip{r) = otherwise, where <I>(r) = 4m [((T/r)-*-^ — (cr/r)^] 
is the genuine LJ potential. The particles were placed 
in a basic cubic box of volume V = L^, and the mod- 
ification of $(r) with fc = L/2 « 3.36(7 as well as 
the periodic boundary conditions have been used to ex- 
clude the finite-size effects. The simulations were per- 
formed in a microcanonical ensemble at a reduced den- 
sity of n* = ^(7^ = 0.845 and a reduced temperature of 
T* = k^T/u = 1.7. All runs of the length in I = 10 000 
time steps each were started from an identical well equi- 
librated initial configuration p(0). The precision of the 
algorithms was measured in terms of the relative to- 
tal energy fluctuations £ = {{E — {E)Y) /\{E)\, where 
E ^ \ Y.f=i ■m'Vi'^/2 + i E^j5^,) Virt]) and ( ) denotes 
the microcanonical averaging. Note that if the equations 
of motion could be solved exactly, the above fluctuations 
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should vanish, because in microcanonical ensembles the 
total energy is an integral of motion, E{t) = E{0). So 
that during approximate MD integrations, smaller values 
of £ will correspond to a better precision in evaluation of 
phase trajectories. 



FIG. 1. The total energy fluctuations obtained in the sim- 
ulations for different values of free parameter ^ at three fixed 
time steps, h* = 0.00125, 0.0025, and 0.005 using the VEFRL 
(subset (a)) and PEFRL (subset (b)) integration (Eqs. (21) 
and (22, respectively). The simulation results are presented 
by circles connected by the solid curves. The function 7(^) 
(see Eq. (13)) is plotted by the dashed curves. 

The equations of motion were first integrated with the 
help of extended decomposition scheme (9) in which the 
parameters ^, A, and x, being constant within each run, 
varied from one run to another. For convenience of pre- 
sentation, we have chosen ^ (instead of A as earlier) to be 
a free parameter, so that the two other quantities A(^) 
and xiS.) should be treated as depending on ^ accord- 
ing to constraints (16). The total energy fluctuations ob- 
tained in such simulations at the end of the runs for three 
(fixed within each run) undimensional time steps, namely 
h* = h{u/ma^y/^ = 0.00125, 0.0025, and 0.005, are 
shown in Fig. 1 as functions of parameter The subsets 
(a) and (b) of this figure correspond to VEFRL (Eq. (21)) 
and PEFRL (Eq. (22)) versions, respectively. As can be 
seen, all the dependencies £{£,,h) have the global mini- 
mum which locates at the same point, £, « 0.164 for VE- 
FRL or ^ « 0.179 for PEFRL, independently on the size 
h of the time step. This point coincides completely with 
the minimum given by Eq. (19) or (20) for the function 
7(^) = 7(^, A(^), xiO) (see Eq. (13), in which we should 
put 74 = and 76 = for VEFRL or 71 = and 73 = 
for PEFRL) which is also included in Fig. 1 (dashed 
curves in the subsets). Moreover, the fluctuations £{^, h) 
appear to be proportional to the norm F = 7/1^ (Eq. (14)) 
of global errors, and the coefficient of this proportional- 
ity almost does not depend on ^ and h. In addition, at 
each time step the energy fluctuations decrease at the 
minimum in 40-60 times (in dependence on the versions 
considered) with respect to those at ^ = 0, that is in 
agreement with our predicted values 7FR/7mkf^^ ~ 43 
and7FR/7^f™^«62. 

The result for the total energy fluctuations as func- 
tions of the length I = t/h oi the simulations correspond- 
ing to the optimized VEFRL (Eqs. (19) and (21))) and 



PEFRL (Eqs. (20) and (22)) algorithms is presented in 
subsets (a), (b), (c), and (d) of Fig. 2 for the time steps 
h* = 0.00125, 0.0025, 0.005, and 0.01, respectively. For 
the purpose of comparison, the same dependencies cor- 
responding to original VFR and PER integrators (see 
Eq. (8)) are also shown there. Note that for the original 
integrators, the time step for each subset was chosen to 
be always in factor 4/3 smaller than that of the optimized 
versions, namely, h* = 0.0009375, 0.001875, 0.00375, and 
0.0075. The last condition is necessary to provide the 
same number of total force recalculations during the same 
observation time t ^ h within both the approaches. Note 
also that within these approaches, the hfth- and higher- 
orders truncation uncertainties become too big at step 
sizes h* > 0.01. In particular, then the ratio of the total 
energy fluctuations to the fluctuations in potential en- 
ergy (the standard quantity to estimate the precision of 
the calculations) exceeds a few per cent. Such large step 
sizes cannot be used in precise MD simulations and, for 
this reason, are not considered here. 



FIG. 2. The total energy fluctuations as functions of the 
length of the simulations performed using the optimized VE- 
FRL and PEFRL algorithms, as well as the original VFR and 
PFR integrators. The results corresponding to different val- 
ues of the time step, namely, h* = 0.00125 and 0.0009375, 
0.0025 and 0.001875, 0.005 and 0.00375, as weU as 0.01 and 
0.0075 are presented in subsets (a), (b), (c), and (d), respec- 
tively. 

As can be seen from Fig. 2, both the original (VFR 
and PFR) and optimized (VEFRL and PEFRL) algo- 
rithms exhibit very good stability properties, that is ex- 
plained by the symplecticity and time reversibility of the 
produced solutions. No systematic deviations in the to- 
tal energy fluctuations can be observed for all the in- 
tegrators. Instead, in each of the cases the amplitude 
of the fluctuations quickly tends to its own value which 
does not increase with further increasing the length of 
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the simulations. However, this value is significantly 
larger within the original VFR and PFR integration. 
At the same time, using the optimized VEFRL and 
PEFRL algorithms allows us to decrease the unphys- 
ical energy fluctuations approximately in factor 10-15 
(for VEFRL) or 20-25 (for PEFRL), despite the larger 
time steps. This is in excellent accord with our theo- 
retical predictions r"^^^^{Ah/3)/rFR{h) w 0.073 and 
rmf„™^(4V3)/rFR(/i) « 0.051, obtained in section IIL 
It is worth pointing out also that the PEFRL algorithm 
is slightly better in energy conservation than its VEFRL 
counterpart (whereas the VFR integrator is better with 
respect to its PFR version) . 

Consider now the results for composition integration. 
The total energy fluctuations obtained in this case at the 
end of the simulation runs for the same set (as in Fig. 1) of 
undimensional time steps, h* = hiu/ma'^Y/'^ — 0.00125, 
0.0025, and 0.005, are shown in Fig. 3 as depending on 
parameters ^ (subsets (a) and (c)) and x ((b) and (d)). 
The subsets (a)-(b) and (c)-(d) relate to the VESL and 
PESL versions, respectively, of composition scheme (25). 
Mention that the time coefficients ^(x) and A(x) for this 
scheme should satisfy equality (31) in which the quantity 
X = 1 — 2(^ + A) is treated as a free parameter. So that 
only one x-dependency is enough, in principle, to present 
the results. But, in view of somewhat complicated struc- 
ture of function £{x), we have included the ^-dependency 
as well for more clarity. Note also that according to con- 
straint (31), the parameter x cannot accept values from 
the interval ]1 — 4'i9 « —0.658, 1], because otherwise this 
will lead to imaginary solutions for ^ and A. 



FIG. 3. The total energy fluctuations obtained in the sim- 
ulations for different values of parameters ^ and x three 
fixed time steps, h* = 0.00125, 0.0025, and 0.005 using the 
VESL (subsets (a) and (b)) and PESL ((c) and (d)) versions 
of composition integration (Eq. (25)). The simulation results 
are presented by circles connected by the solid curves. The 
function 7 (see Eq. (29)) are plotted by the dashed curves. 



As can easily be seen from Fig. (3), all the dependen- 
cies have up three minima and their locations on ^- and 
X- axes practically do not depend on the size h of the 
time step. Moreover, the fluctuations £(S,, h) again ap- 
pear to be proportional to the norm F = 7/1^ of global 
errors, where the function 7 is now deflned by Eq. (29) 
and plotted in the flgure as well. Among the three min- 
ima, two of them have nearly the same depth, and for 
VESL integration the global minimum of £ is achieved 
at ^ « 0.323 or X ~ —.726, i.e., from the left of the in- 
terval ]1 — 4i?, 1] for X- This coincides with the predicted 
values given by Eq. (35) corresponding to the global min- 
imum of 7. In the case of PESL integration the pattern is 
somewhat different. At small enough values of the time 
step, namely at h < 0.0025, the global minimum for the 
functions £ and 7 is identified from the right of the above 
X-internal, namely, at ^ w 1.45 or x ~ 2.4. But with fur- 
ther increasing h, the position of this minimum for £ 
begin to shift with respect to its initial value. Such a be- 
havior can be explained by the influence of higher-order 
0(/i^) truncation errors, which appears to be significant 
in this case due to the largeness of absolute values of 
time coefficients ^, A, and x (in particular, then all these 
coefficients are larger in amplitude than unity). This is 
contrary to a left-lying minimum of £ which is achieved 
at the same point ^ « 0.316 or x ~ —0.737 as that of 
function 7 (see Eq. ((34) independently on the size of the 
time steps considered. Then taking into account that 
both the minima have almost the same depth, and the 
fact that the left-lying minimum of £, being local at small 
h, transforms into the global minimum at moderate and 
larger time steps, solutions given by Eq. (34) should be 
considered as optimal. 

As expected, at the minima found the energy fluctu- 
ations decrease approximately in 1.5 times with respect 
to those at = A = I? w 0.41 (and x ~ —0.66) cor- 
responding to the original Suzuki approach. It is inter- 
esting to point out that the constraint ^ = A minimizes 
the function J2q=i M9I = + 2|A| -I- |x| and provides 
a minimum for the maximal value among the quantities 
1^1, |A|, and |x|- Using these last two criteria and the 
additional requirement that \dq\ < 1, Kahan and Li |l^] 
have derived extended composition schemes like (23) up 
to the tenth order in the time step. In view of our more 
precise analysis performed in this paper we see, therefore, 
that the previous criteria for optimization of algorithms 
may have nothing to do with providing the best perfor- 
mance of the computations. Finally, comparing the re- 
sults presented in Figs. 1 and 3 between themselves, it 
can be pointed out that the values for minima of £ corre- 
sponding to the VESL and PESL integration are approx- 
imately in 1.5 times higher than those of the VEFRL and 
PEFRL algorithms. In addition, VESL and PESL inte- 
grators require more (five, instead of four) force evalua- 
tions per time step, so that their efficiency is in factor 1.5 
(5/4)4 _ 3 ^Qj.gg ^^^^ respect to the VEFRL and PEFRL 
schemes. Therefore, as was predicted (see. Eq. (33)), the 
preference should be given to the decomposition VEFRL 
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and PEFRL algorithms rather than to their composition 
counterparts VESL and PESL. 

VI. CONCLUSION 

We have formulated a new approach for improving the 
efficiency of decomposition integration of the equations 
of motion in many-body systems. As a result, the Forest- 
Ruth-like algorithms have been derived and optimized for 
MD as well as celestial mechanics simulations. The algo- 
rithms are explicit, simple in implementation, and pro- 
duce phase trajectories which are exactly symplectic and 
time reversible. Their main advantage over the original 
integrators by Forest and Ruth is the possibility to gener- 
ate much more precise solutions at the same overall com- 
putational efforts. It has been predicted theoretically and 
confirmed in actual MD simulations of a L J fluid that the 
optimized Forest-Ruth-like algorithms allow to decrease 
significantly the unphysical fluctuations of the total en- 
ergy using even greater time steps. The question of how 
to optimize the integration by composing second-order 
schemes has also been studied. As has been shown, the 
Suzuki-like algorithms, obtained in this case, although 
improve the efficiency of the computations with respect 
to the original integrators by Suzuki as well as Forest 
and Ruth, appear to be less efficient than the optimized 
Forest-Ruth-like algorithms. The last algorithms should 
be considered as optimal among all possible decomposi- 
tion integrators of the fourth-order with single splitting 
of the Liouville operator. 

The approach presented can also be adapted for op- 
timization of the integration of motion in more compli- 
cated systems for which splitting of the Liouville operator 
into more than two parts may pay dividends. The exam- 
ples are multi-component fluids and systems with long- 
range interactions, where characteristic time intervals of 
the dynamical processes can differ by many orders from 
each other. In some cases, for instance, for systems with 
orientational degrees of freedom, the additional splitting 
may be necessary to obtain analytically integrable parts. 
These and related problems will be considered in a sepa- 
rate publication. 

Part of this work was supported by the Fonds 
zur Forderung der wissenschaftlichen Forschung under 
Project No. P15247. 

Appendix. Optimization of Forest-Ruth-like 
algorithms in specific cases 

In section III, we have derived the velocity- and 
position- Forest- Ruth-like algorithms and optimized them 
for integration of motion in MD simulations. Now, one 
considers two speciflc cases when the efficiency of decom- 
position scheme (9) can be improved additionally. One 
example is a collection of weakly interacting particles, in 
which the potential part B = eB oi the Liouville opera- 
tor can be treated as a small perturbation with respect 



to its kinetic counterpart A, where e -^1. Another case 
is the integration of motion in the solar system, where 
the interactions between planets are negligible small in 
comparison with the force that acts on them with respect 
to the sun. Then, the Liouville operator can be presented 
as i = A + eB, where now A is the sum of two-bodies Li- 
ouville's commutative operators, each of which describes 
the motion of an isolated body in an external field of the 
central mass, whereas sB corresponds to weak interac- 
tions of the bodies between themselves. Here, as in the 
usual case, the separate exponential operators e-^"^ and 
e^'^ appear to be exactly integrable, so that the decom- 
position scheme will also lead to explicit integration. 

Therefore, assuming that the parameter e is sufficiently 
small and considering the velocity version of extended 
Forest- Ruth- like decomposition (9), we can neglect by 
fifth-order commutators arising in the truncation term 
C5 at the multipliers 72, 73, 74, 75, and 75 (see Eq. (12)). 
The reason is that such commutators are of order or 
higher, whereas the commutator at 71 is proportional 
only to the first power of e. In such a case, it is quite 
natural to reduce the multiplier 71 to zero exactly, rather 
than to find a minimum for the norm (13) of 7 with re- 
spect to all the commutators. Then the problem (15) 
transforms into the following system of three equations 

' a(^, A,x) =0, 

< /3(e,A,x)=0, (Al) 
Ji(^,A,x) = 0, 

where a and (3 are defined according to Eq. (11), and 
the expression for 71 follows immediately after Eq. (12). 
Since there exist up four sets of real solutions to Eq. (Al), 
an additional requirement should be imposed to choose 
the optimal values among them. Obviously, such values 
must minimize the norm of the main term of truncated 
uncertainties, which now is equal to VtF+Ts (note that 
commutators at 72 and 73 are proportional to e^, while 
to the third power of e or higher at 74, 75, and 75). As 
a result, one finds the optimal solutions 

(, = 0.5316386245813512 

A = 0.5437514219173741 (A2) 

X= -0.3086019704406066. 

For the position version (when A ^ B), the commuta- 
tors at 71 and 74 are interchanged, so that it is necessary 
to solve the system a = 0, /? = 0, and 74 = 0. No real 
solutions to this system have been found. In such a sit- 
uation, we could try to find the global minimum for I74I. 
But since it appears to be greater than zero, the position 
version is not recommended to be used for this case and 
the preference should be given to the extended integra- 
tion in velocity form (with time coefficients deflned by 
Eq. (A2)). Note that the main terms of uncertainties for 
this integration are 0{e^h^) and 0{eh'), so that in view 
of the smallness of e, the scheme derived can be related 
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to a quasi sixth-order algorithm. Such an algorithm, con- 
trary to usual sixth-order schemes, will require four, in- 
stead of seven [§,0, force evaluations per time step. 

In the above derivation, we tentatively assumed that 
the order of smallness of e is nearly the same as that 
of the size for the time step h. In systems for which 
the parameter e is extremely small, the strategy of opti- 
mization should be reconstructed. Namely, we can now 
ignore all the terms of the second order with respect to 
e, even that which appears at /3 (see Eq. (10)) in third- 
order uncertainties Oilp?) with respect to h. Then within 
the velocity version, one comes to two equations, a — Q 
and 7i = 0, for three variables, ^, A, and Thus, the 
order of the decomposition scheme can be increased addi- 
tionally by canceling an extra-order term OielJ) (which 
is assumed to be much more significant than the terms 
0{e'^h^) and 0{e^h^)). An exphcit expression for 0{£h7) 
(it is denoted simply as 0{h7) in decomposition formula 
(9)) has the form 

Oieh') = (4^^ [A, [A, [A, [A, [A, [A, B]]]]]] + 0{e^))h^ , 
where 
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So that the desired system of equations is 



A,x) 
7i(C, A,x) 

4''^(e,A,x) 



(A3) 



It can be solved analytically, and the solutions are 



^ 20' 



2 V 7 



X 



— , (A4) 
180 ' ^ ' 



where the preference should be given for sign be- 
cause this leads to a smaller value (« 0.0036) for (and 
for an multiplier at the ninth-order commutator) which 
forms the main term 0{£^h^) (and 0{£h^)) of local trun- 
cation errors. 

In the case of position version (A ^ S), the seventh- 
order (with respect to h) truncation term behaves as 

0{eh^) = (4^)[i?, [S, [i?, [B, [i?, [A, B]]]]]] + 0{e^))h' , 



with 
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and equations (A3) transforms into 

/?(e,A,x) = 
74(e,A,x) = 

4''^(^.A,x) = 



(A5) 



System (A5) has eight solutions and all of them are 
real. The optimal solution, which minimizes (to the value 
w 0.0034) the multiplier \a\ (it forms now the main term 
0{e'^h^) of local uncertainties), is 
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30 



(A6) 
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For both velocity (A4) and position (A6) versions, the 
main terms of uncertainties are 0{e'^h^) and 0{eh^). So 
that, taking into account the smallness of £, these ver- 
sions present, in fact, quasi eight-order algorithms, which 
contrary to usual eight-order schemes, will require again 
only four, instead of up fifteen force evaluations per 
time step. 
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